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| Section 1 

» INTRODUCTION AND SUMMARY 

1 This contract was modified on 13 April 1972 to include an extended work 

statement and period of performance through 13 January 1973. This final 
] report documents results achieved during this extended effort. The additional 

effort consisted of two tasks, each of which has been separately documented in 
] Lockheed technical reports that are appended to this report. The two tasks 

are described as: (1) "Backflow of Outgas Contamination onto Orbiting Space- 

] craft as a Result of Intermolecular Collisions," and (2) "Particle Cloud Buildup 

* Due to Waste Dumping from Orbiting Spacecraft." Both of these tasks deal with 

| aspects of spacecraft contamination in orbit. 

? Spacecraft contamination is a matter of current and vital interest in the 

j design of thermal control surfaces, solar panels and experiment and instru- 

mentation packages for use aboard earth-orbiting vehicles. Contamination 
| sources include liquid waste dumps, reaction motor firings and material out- 

gassing. The two tasks performed under this contract modification provide 
analytical techniques with which quantitative estimates can be made of certain 
specific contamination threats. 

! Some of the outgas products flowing away from a spacecraft will be de- 

flected back onto the spacecraft as a result of intermolecular collisions. This 
j is particularly true with respect to collisions with ambient atmospheric mole- 

cules. Under Task 1, analytical techniques were developed to predict the amount 
of mass that will return. For the Skylab spacecraft, it was estimated that on 
the order of 1 or 2% of the outgas products will return to the spacecraft as a 
result of intermolecular collisions. 

Mass dumped overboard frcm liquid waste dumps will leave the space- 
craft and go into orbit about the central body (Earth, Moon, etc.) in much the 
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same fashion as the spacecraft. The effect of continuing waste dumps is to 
create a particle cloud in orbit that continuously expands along the spacecraft 
orbit. Some of this mass will return to the spacecraft at points where the 
particle orbits intersect the spacecraft orbit. It was estimated for the Skylab 
spacecraft, neglecting aerodynamic drag, that a particle cloud will extend, on 
both sides of the spacecraft, about 36 km in the radial direction, 9 km in the 
out-of-plane direction and will expand along the spacecraft orbit at the rate of 
168 km for each spacecraft orbit. Aerodynamic drag will cause the particle 
cloud to spiral inward toward the central body. Again neglecting aerodynamic 
drag, the estimated amount of return mass flow was found to be about of the 
same order as a typical spacecraft material outgas rate. When aerodynamic 
drag was considered, it was found that particle collisions with the spacecraft 
were impossible, and the return mass flow rate was, therefore, zero. The 
predicted light scattering properties of the optical cloud were found to be 
nearly identical with those predicted for an isotropic source, with some in- 
tensification in directions along the spacecraft orbit. 
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Section 2 

DISCUSSION AND CONCLUDING REMARKS 


The results of the Task 1 effort, "Backflow of Outgas Contamination 
onto Orbiting Spacecraft as a Result of Inter molecular Collisions," is docu- 
mented in Ref. 1 and is included in this final report as Appendix A. As pointed 
out in Appendix A, the treatment of backscatte ring due to collisions between 
outflow molecules and atmospheric molecules is considered to be fairly 
rigorous. The treatment involving collisions between individual outgas mole- 
cules, however, is based on somewhat extreme simplifications. Although the 
backflow due to self-scattering was found to be small for the Skylab spacecrait, 
future applications might make worthwhile an effort to develop a more ngorous 

treatment. 


The second task, "Particle Cloud Buildup Due to Waste Dumping from 
Orbiting Spacecraft," is documented in Ref. 2 and is included in this final report 
as Appendix B. There are a number of areas where this study might be expanded. 
In the calculations of the particle cloud buildup, it was assumed that the expulsion 
rate was continuous and uniformly distributed in all directions. It was furtner 
assumed that the particle size and velocitv distributions were uniform and constant 
with respect to time. A worthwhile future effort would be to remove some of these 
restrictions and also to consider the long-term interaction and thermal equilibrium 
of these particles with the ambient atmosphere. 
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FOREWORD 

This interim report documents the results accomplished 
during performance of Task A of Contract NAS8-26554, "Mass 
Transport Contamination Study," Modification 1. The work was 
performed by personnel in the Aeromechanics Department of 
the Lockheed -Huntsville Research & Engineering Center under 
Contract to NASA-Marshall Space Flight Center. The MSFC 
technical monitor for the contract is E.E. Klingman, S&E- 
SSL-PM. 
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Section 1 

INTRODUCTION AND SUMMARY 

Orbiting space vehicles lose mass to the ambient environment as a re- 
sult of outgassing cf spacecraft materials, dumping of waste products and 
firing of attitude control reaction motors. Some of the lost mass will return 
to the vehicle, where even very small amounts may significantly cloud space- 
craft windows, contaminate experiment packages or otherwise interfere with 
spacecraft operations. 

In the highly rarefied environment at orbital altitudes, backscattering 
due to intermolecular collisions will be a primary mode by which outgas 
products will be deflected back onto the vehicle. The intermolecular collisions 
will be of two distinct types. First, they will occur between the outflowing 
molecules and the ambient atmospheric molecules. The stream of ambient 
atmospheric molecules will be highly energetic, corresponding to the orbital 
velocity of the spacecraft. Intermolecular collisions of this type will conse- 
quently tend to strongly deflect the outflowing molecules in the direction of the 
ambient freestream. Second, intermolecular collisions will occur between 
individual outflowing molecules because o. differences in velocity (including 
direction) of the individual molecules. 

To provide '.n indication of the amount of backflow to be expected, a 
theoretical analysis was performed based on the following simplifying assump- 
tions. The spacecraft was assumed to be spherical in shape with the mass 
flow emitting uniformly over the spherical surface at a constant rate and in a 
D’Uambertian spatial distribution. The outflow gases were assumed to be 
neutrally charged and of a single species with a molecular weight character 
istic of a composite cf the actual species involved in the mass flow. 
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The theoretical analysis showed that, tor outgassing only, less than 1.5% 
of the outgas products will return to the Skylab spacecraft as a result of inter- 
molecular collisions. When the total mass flow from the spacecraft, including 
waste dumps, reaction control motor firings, etc., was considered, it was esti 
mated that about 30% will return to the spacecraft. The latter result is based 
on certain rather extreme simplifying assumptions. 
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Section 2 

METHOD OF APPROACH AND DERIVATION OF EQUATIONS 
2.1 METHOD OF APPROACH 

The highly rarefied environment under consideration is such that the 
flow may be considered to be "nearly free molecular." In this flow regime, 
the mean free -path between intermolecular collisions is large compared to 
typical spacecraft dimensions, yet small enough to cause significant departure 
from free molecular flow. The fact that the flow regime is indeed nearly free 
molecular will be shown later for the Skylab spacecraft. 

In this flow regime, the "first-collision" model approach of Baker and 
Charwat (Ref. 1) is appropriate. Each molecule is considered to undergo a 
maximum of one intermolecular collision prior to colliding with a surface, or 
between two successive collisions with a surface. In the present applications, 
this allows two classes of collisions to be considered: (1) between an outflowing 
molecule and an ambient atmospheric molecule, neither of which has undergone 
previous collisions; and (2) between two outflowing molecules which have not 
experienced prior collisions. Collisions between ambient atmospheric mole- 
cules are neglected altogether since the frequency of such collisions is ex- 
tremely small compared to that of the other two classes of collisions. Actually, 
Baker and Charwat considered only class (1) collisions. Since it is difficult to 
show that class (2) collisions are in fact negligible (Ref. 2), they will be con- 
sidered in this analysis by using certain simplifying assumptions. 

The first-collision approach considers the molecules to be smooth, hard 
elastic spheres of some finite diameter. In the present analysis, the diameter 
is made flexible depending on the relative velocity Detween collisions to account 
for realistic temperature variations in viscosity. 
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2.2 DERIVATION OF EQUATIONS 

It is necessary first to determine the collision rate per unit volume that 
occurs in the region surrounding the sphere. Next, a determination is made 
of the number of colliding outflow molecules in an element of volume which are 
deflected back in the direction of a given point on the surface. The total back- 
flow flux at the given point is then determined by integrating over the half space 
outward from the point on the surface. This is described separately for the two 
classes of collisions in the following paragraphs. 


2.2.1 Scattering by Ambient Atmospheric Molecules 


Consider the geometrical representation shown in Fig. 1. The point on 
the sphere surface at which the backflow flux is to be calculated is designated 
Q, while a point in the field surrounding the sphere which is a source of back- 
scattered molecules is designated P. For convenience Q is considered to be 
in the T, plane, where i, j and k are unit vectors in the direction of corre- 
sponding rectangular coordinates with i pointing in the negative direction of 
the ambient flow. The angle is the latitudinal displacement from the i 
axis. In the primed coordinate system, i' is in the direction of the radius 
passing through Q, j* is in the i, j plane and k is in the k direction. <j> and 
0 / are polar coordinates about the i axis. Finally, r is the distance from che 
sphere center out to P, while r' is the distance from Q to P. 


The collision rate per unit volume, n, in the field surrounding the sphere 
is given by: 


n 


n 


e 


n 

oc 




a 

e -oo 


( 1 ) 


where n is the density field of outflow molecules, is the density of the 

ambient atmosphere, is the ambient flow velocity, v^ is the velocity u* 

the outflow molecules and a is the cross section for collisions oetwcen 

e -oo 

the outflow molecules and the ambient atmospheric molecules. 
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The density field n of outflowing molecules at P is determined by 


f Sv dw 

n e = / ~ ~ 

a w 

f2 


( 2 ) 


where q is the outflow flux at the sphere surface, v is the mean velocity 
w w 

of the outflow molecules and the integration is over the solid angle w sub- 
tended by the sphere at P. We assume here that the variation in the density 
field is due entirely to the diverging flow away from the sphere, and that 
attenuation due to collisions may be neglected. This assumption is justified 
by the mean free-path being much larger than the sphere radius, as is shown 
later for the Skylab spacecraft. 


Carrying out the integration in Eq. (2) yields 

2 q S 2 IT 

w r - \r - R 


w 


Referring again to Eq. (1), n is seen to be given by 


(3) 


n = 


2 q n 
w oo 


V - v O 
oo e e -oo 




W 


(4) 


Next, the directional distribution must be determined of scattering of 
those outflow molecules undergoing collisions in an element of volume ai P. 

To do this, the velocity of the outflow molecules is assumed to be negligible 
compared to the velocity of the ambient atmospheric molecules. The outflow 
molecules may then be considered stationary with respect to the sphere. Con- 
side! the hard sphere collision shown schematically in Fig. 2. For all collisions 
in which the line of centers is in the solid angle element dw, the stationary mole- 
cule will be deflected into dw on the opposite side of the stationary molecule. 

The probability dP^ of the line of centers lying in dw is given by the amount 
of sphere surface contained in dw projected in the direction of flow ratioed to 
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- Diagram for Scattered Flux Due to Collisions 
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the sphere cross section: 

dP 


w 


z 

a du cosa j 

- — ^ cosa = du) 

2 7 r 

7T a 


(5) 


Now, referring to Fig. 3, the flux d q fa is found of scattered outflow 
molecules from collisions in volume element dV at P and crossing surface 
element dA at Q. Surface element dA is oriented at angle fi to the line join- 
ing P and Q, and Q is at a distance r' from P. 


dc U 


. ,,, cosa j. . 
n dv du 

7T 

r' 2 dw 


^ h coso cos (3 j,, 

cos|3 = 2 dV 

7T r' 


( 6 ) 


Referring again to Fig. 1, the flux dq b of backscattered outflow mole- 
cules at Q is 


2 q n V a 
w oo oo e -oo 

C cosa cos <t> 1 

- >1 

r 2 - R 2 \ 

q b * ” w J 

1 r ,2 1 

V ' / 


dV 


(7) 


where the volumetric integral is carried out over the portion of the half space 
outward from dA at Q which is also forward of Q in the negative direction of 
the external flow (recall from Fig. 2 that molecules are scattered only in the 
forward direction). 

Before Eq. (7) is evaluated for arbitrary <f> , the integration is performed 
at <t> = 0 and at <t> = jr/2. 


i, 4 n V a / 7 

b >° _ °° oo e-oo I cog 2 

q VI 

H w w 


tt/ 2 oo 

j cos 2 *' .in * f r or'd*' 


( 8 ) 


The radial distance r under the r' integral sign is a function of r' and V 
as follows 


r = (R 2 + 2 cos <f> r' R + r' 2 ) 
2-6 


1/2 


( 9 ) 
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The difficulty of performing the integration is lessened greatly by approxi- 
mating the expression under the r' integral sign by the following simpler 
expression: 


r - Vr 2 - R 2 ~ 


~ J_ (*} R^\ 

= 2 \r 2 r 6 / 


( 10 ) 


The close agreement between the two expressions is shown in Fig. 4. 


Equation (8) now becomes 


2 n V a R 
b, o _ oo oo e -oo 


'w 


w 


it / 2 

J* cos 2<f>' sin$ 


Jf — *1——^ f - 

\J 1+2 co8(f> r + r J (\ + 


dr' 


~ (\ + 2 CO80' *?' +7'^! 


3-]d0 (11) 


where r' = r'/R. The group of constant factors preceding the integral signs 
may be rearranged and expressed in terms of certain dimensionless param- 
eters widely used in rarefied gas dynamics: 


2 n V a R 
oo oo e -oo 

V 

w 


S, a 


= Jz 

jf 2 Kn 


e -oo 


Kn a 


( 12 ) 


oo 


* 


where S, is the ambient velocity ratioed to the most probable thermal velocity 
b 

of the outflow molecules at the surface temperature, Kn is the ambient mean 
free -path ratioed to the sphere radius, and cr ° 88 section for 

collisions between outflow molecules and ambient atmospheric molecules 
ratioed to the cross section for collisions between individual ambient atmos- 
pheric molecules. 
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3 4 

r/R 


5 6 7 8 9 10 


2 2 6 ^ 

- Approximate Expression R /r + R /r Compared to Exact 


Expression 2(r - \r^-~R^)/r 
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The integrals over r' in Eq. (11) are given as follows: 

o© 


/ 


dr' 


B _X- 

1+2 cos <j> 7 , +'r'^ Sin ^ 


(13) 


00 


A 


dr 


3$' cos^' 3 cos<ft 


i + 2 cos <t> r ' + r ' 


,2\ 3 8 sin 5 <t>' 4 sin^V 8 sin**^' 


. 2 . 


. 4 


(14) 


Equation (11) now becomes: 

Vf 


1 b, 0 
q w 


S, o 
b e -oo 

Kn a 


E /2 


J* <p cos^^d^ 


/ / 2 


./ 2 . » 

9 c .°g--£ d^' 


• 4 a' 
sin 9 


jt/2 


ff /2 


i / -!/ 

Vf 


3 i/ 

COS 9 jjl' 

_ V * 

sin (p 


37T 2 - 8 a e-oo 


Kn o 


(15) 


For <j> = ir/2, making use of the results of Eqs, (10) and (12) through (14), 
Eq. (7) becomes: 


q b, 7 T / 2 _ 1 fb CT e -oo 


w 


V2? 


Kn a 


oo 


jr/2 7r/2 

cos9 d0 sin^'cos 0' 


» 3 1 cob 6 ' 3 cos<ft' 

, ^ + 8 . 4., “ 4 sin^' 8 .3., 

j sin ^ T sin 0 


df = 


ft 


S, ° 
b e -oo 


32 12 Kn a 


( 16 ) 


00 


or 


Lt W-j = = 0.3471 

q b,0 2(3 jr - 8) 


(17) 
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For general^, the integration of Eq. (7) becomes much more tedius. The 
angle a is expressed in terms of ^ , <f>' and 0 as follows: 


a = cos' 1 (cos <p cos - sin<P sin$ cos9 ) 


Recall that the volumetric integral in Eq. (7) is carried out over the portion of 
the half space forward of Q in the negative direction of the external flow. This 
is equivalent to requiring that 0 <. a < jr/2. For 0 < $< jr/2, Eq. (7) becomes: 


1 S u 

1 b e -ao 

=— Kn a 


!-♦ / 


J* (cos <(> cos<f> - sin <f> sin<f>' cosQ') cos<f>' 


„ (a + 1 £- . 1 . 2 _COs£\ d9 > 

x i + 8 . 4 4 sin^' 8 . 3 .>) aw v 

V sin to r sin <p / 


sin q> 
it/Z rr 


f f, , 


(cos <f> costf)' - sin<£ sin<// cosQ*) cos^ 


tan<£ tan^> ; 


x /a' + 2 . l££2l' . | d6 ' dA 

v s .4., 4 sin^ ' 8 . B 3j v 

\ sin 0 sin 0 ' \ 


Performing the integrations over 0 reduces Eq, (^9) to: 


i S, a 

J b e -oo 

x- Kn a 


x in-cos^ 


V L’ +| - I— 1 d*' 

* T 8 . 4 ,» 4 sin9 8 . 3 /I 

' Sin 0 sin (ft/ 


8 • 3 .» 

sin f 


cos^ 


C08 ‘WW C08 ‘ 


IV .3 1 CO 9^ C03^' \ d v 

* r 8^T7‘ 4if^ .3.')“* 

\ sin 9 sin 9 ' 
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+ cos 


♦ t 1 cV - H$'- *■?£)<♦' 

J 1 ' sxn 9 sin 9 


jr 
2 ‘ 


The first integral in Eq. '20) is identical to the integral in Eq. (15): 


jr/ 2 

/ 


2 j 

COS 9 


/ /a' , 2 . i £°2l . 2-£°2£\ <39 = 

r + 8 . 4 4 sin^ 8 . 3 ,'/ 

' sin 9 Bin 9 


3r 2 - 8 
48 


( 20 ) 


( 21 ) 


After some algebraic manipulation is performed, and expressing in terms of 
the flux at <t> = 0, Eq. (20) then becomes: 



To allow integration in closed form, the first factor under the integral 
sign is approximated as follows: 



tan 2 (| - 9 ) 
tan 2 ^ 


tan(-|-?) tan(|-- 9 ) 

tan? COS tan? 


= 1 


tan(?- - ) 

I tan? 



tan 2 (f-?') 

2T 

tan 9 


(2 j) 


The approximation was obtained by fitting the quadratic equation to mate’ 
end points at x = tan(| - ? )/tan? = 0 and 1 and the first derivative at x - 0 
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Very close agreement between the exact and approximate expressions if- shown 
in the comparison in Fig. 5. 

Equation (22) may now be expressed as follows: 


jt/2 



The first integral is seen to be identical to the integral in Eq. (16) for 
4 = */2, while the other two integrals drop out due to the cos 4 factors. Also 
all three integrals drop out at 4 = 0 due to the identity of the limits of inte- 
grations, leaving the first term equal to unity. At 4 = 0 and jt/ 2, therefore, 
Eq. (24) reduces to the previous results for those particular points. 

Carrying out the integrations in Eq. (24) yields: 

q h / , 48 

— — = Ccos9 + 5 

q b.o \ 1(3/ -8) 

- i sin24 +| (|tan 2 4 - — V" * tan *j] 

8 8 6 cos 4 ^ 



sin 


' 4+4 
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-|co 9 * (§-*) sin24> -* 2 - \ain Z <t> + £ tan 2 * (<§ - <M tan* - l)j 


2 , r 


+ (f ' ^70 > 37 ( 3( l - tan ^ - 3 tan ♦ + tan * 


+ 2 sin2* - 5* 


>j - - *) lncos* ~ 2 ~ lncosx dxi 


There remains now only the task of evaluating the integral of the lncosx. 
To accomplish this, lncosx is approximated as follows: 


~ i / 1 ~ X \ + i x 7;2 ~ - 4 x 2 
'.rxofx = i n l _ ) + _ x " 2 

\ f / * Z* 


This expression was obtained by matching the end point values and the first and 
second derivatives at x = 0, and requiring the limit of the ratio of the two ex- 
pressions to approach unity as x — tt/ 2. Good agreement is shown in the 


comparison 


of the two expressions in Fig. 6. Performing the integral yields: 


lncosx dx = • | id - ~ V) Ml +f <> + 


l aZ r - t 


* 


Combining this result with Eq . (25) yields th< 


tinai res 


' jr(3 V - 8) ( l. 


- i sin2* + | (“ tan 2 * - — K- - tan * ) 
& 82 cos^ 


f— ^ Z 

^ cosp | it* - s tn2^> - C - ^ ain 0 
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•0 


+ \ tan 2 $ ((§ - 

> 

+ 2 sin2* - 5* 


) tan* - l)j + (f - 1 ) 3 ^ (3(f “ ♦) tan ^ “ 3 tan ^ + tan< * 

- (| - <f>) lncos <t> -j(j -4) sin 2 <f> +f |(1- % 4>) ln(l -|tf) + “ *) 




<t> 3 


> 


(28) 


For 7r/2 < (j> < n, Eq. (7) becomes 



1 

V2ff 


a 

e -oo 


Kn a 

00 


-cosy 


7t/2 


_ j i 2 

cos — TTrrrrXf cos 


ff/2-</ / 

7r/2 


cam// tan^' 


,/ L' , 3 <t> I cos^ _ 2 COS^ \ d ( 

* r + 8~4-, • 4iiUf 8 Q ,jJ 9 ' 

\ sin 6 sin 0' 


+ cos 


ip I *an 2 0 tan<*>' - 1 cos V /<£ +| ~^T,' 4 flmj/ ' 8 ~3^') d * 

J . , J \ sin $ sin 


Jr/2-l/' 


(29) 


where ip = ir - <j> . Comparing Eq. (29) to Eq. (20) reveals the two equations to 
be identical except for the absence in Eq. (29) of the first integral in Ec. (20), 
and the substitution of n - <£ for 4> . The results obtained for 0 < <f> < n/2 can 
be used, therefore, to obtain resuits for ;r/ 2 < <t> < ir : 


- ^hJLlSt . cos ^0) 
^b, o ^b, o 

Note that the results for > tt/ 2 ignore the absence of scattering, ’.a 
the "shadowed’ 1 volume behind the sphere. This is Oi no real consequence, 
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however, since, as is shown in the following paragraph, the calculated back- 
scatter flux rapidly decreases as <f> increases beyond 7r/2. 


The variation in backscatter flux with <f > , as determined from Eqs. (17), 
(28) and (30), is given in tabular form in Table 1, and graphically in Fig. 7. 


Table 1 

VARIATION OF BACKSCATTER FLUX WITh ANGULAR DISPLACEMENT 
ON SPHERE (INTERACTION WITH AMBIENT ATMOSPHERE ONLY) 
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2 2.3 Scattering Due to Collisions Between Outflow Molecules 

It will be assumed for the moment that the contribution to the backscatter 
flux due to collisions between individual outflow molecules is small compared 
to that due to collisions with ambient atmospheric molecules. This assumption 
will be shown later to be essentially true when outgassing alone is considered. 
The truth of this assumption justifies the use of certain simplifying assumptions 
in the derivation of the backflow due to self -scattering which would otherwise 
be considered somewhat extreme. The assumptions will produce conservative 
results in that the calculated values will be greater than would be the case with- 
out the assumptions. In addition to the usual fir st -collision assumptions, it is 
assumed that the outflow molecules contained in an element of volume are com- 
posed of two groups, with the two groups traveling in opposite directions at 
their mean thermal velocity. The collision rate per unit volume is determined 
from the collisions between these two groups. In addition, it is assumed that 
the scattered molecules resulting from these collisions flow outward from the 
volume element with a spherically symmetrical directional distribution. 


Based on the above assumptions, the collision rate h per unit volume is 
given by: 


n 


n n 

- — — • — - • (2 v ) o 

7 7 v w' e-e 


(3 A ) 


where a is the cross section for collisions 

e -e 

Making use of Eqs. (3) and (10) yields: 


between outflowing molecules 


n 



o 

e -e 


v 

w 




(32) 


Now, referring again to Fig. 3, the flux d q b of backscattered molecules f.orr. 
collisions in volume element dV at P and crossing surface element dA a. C s 


d q 


b 


n cosb 

4 ' 2 

4 7T r 


(33) 
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Using Eq. (32) and integrating over the entire half-space outward from dA at 
a point Q on the sphere surface (refer again to Fig. 1) yields: 


q, q a R 
w e-e 


w o 


. / f dr' 

sin<p cos<p I — 2 

O ^1+2 cos <f> r' + r ,2 ^ 


. r dr' . r dr 

2 / 4 * + / 1 

o (l +2cos^r'+ r' 2 ) 0 (l+2cos/ r'+r' ) I 


d (ft (34) 


The grouping of parameters preceding the integral sign in Eq. (34) may be 
considered an inverse Knudsen number characteristic of self-scattering . 


q a R ! 

w e ~ e - = Kn' 1 
— e-e 


The integrals over r' in Eq. (34) are as follows: 


1 + 2 cos (ft r' + 7' 2 j 


J_ (ft' ]_ ccs<t> 

2 ~ 2 . 3 2 . 2 

sin (ft sm (ft 


dr' 5 <j[ 5 cosl 5 cosift J_ cos£ 

73 ~,2\ 4 2 16 sin 7 *' ' 24 .inV 16 sinV 6 ^ 2 i’ 
1+2 cos<f> r + r I 


+ 2 cos^'r'+r' Z | 


63 co s (ft 63 cost 


256 . 11 V- . 8 ' 256 . 10/ 

sin (ft sin 0 sin <p 


21 co&6' 9 vQsC 1 cosi 


160 . 6 ‘ 80 . 4 • ,0 .2/ 

sin ^ sin (ft sin (ft 
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Using Eqs. (36) through (38), Eq. (34) now becomes 


= Kn 


'w 


-1 

e -e 


7r/2 


7i’/2 




COS 0 ■, j.' 

— 5 / d< p 

sin 0 

rr/2 


f/ 2 , j/ 2 2 / // 2 

. 5 / / cos(i , a / 5 / cos 0 j/ 5 / 

+ 1 / ft** 'Til ~ v d ^-s / 

V sin •/ sin <p y 

//2 7 // 2 , l /Z 2 / 

1 / COS i 1 / . 63 / x r CO t(fr J J ' 63 / COS ,/ 

-J *25bJ -7j d * 

00 0 

f /2 ? , //2 2 , j / 2 2 / 

63 / cos> _21_ / cos 0 d v 9. / cos .4 d 0' 

' 256 / 9/ d<? ‘ T60 / • 5 ,/ d0 80 / . 3 ./ d9 

CDO J sin 7 0 «/ sin 0 r' sin 0 


ir/2 


-w/ ^' d *i 


(39) 


Carrying out all of the integrals in Eq. (39) would be extremely laborious 
and, considering the assumptions involved in the derivation of this equation, 
rather unprofitable. For this reason, only the first term within the brackets 
in Eq. (34), and, hence, the first two terms in Eq. (39), will be retained in the 
following development. Carrying out the integration of these terms yields 

%. = iii. Kn’ 1 (40) 

q 4 e-e 


* 


I 

1 
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Section 3 

APPLICATION TO FKYLAB SPACECRAFT 

A typical configuration of the Skylab spacecraft, with the Apollo Tele 
scope Mount (ATM), is shown in Fig. 8. Ai equivalent size sphere would be 
roughly 15 to 20 m in diameter. For purposes of this analysis, an equivalent 
diameter of 20 m (or a radius of 10 m) is assumed. 

The operating altitude of the Skylab spacecraft is 435 kilometers. At 
this altitude, the ambient mean free path is 1.38 x 10 4 m (Ref. 3), and the cir- 
cular orbit velocity is 7640 m/sec. The spacecraft surface temperature will 
vary widely, ranging between about 235 and 380°K, depending upon position in 
orbit and other parameters. Assuming a surface temperature of 300°K and a 
molecular weight of 24, the most probable thermal velocity of the outflowing 
molecules is 456 m/sec. Based on the above data, the ambient Knudsen number, 
Kn, is 691 and the speed ratio S b is 16.7. 

The collision cross section <T e-fl0 will vary depending on the relative 
velocity between the colliding particles. This variance is reflected in the 
temperature variation of viscosity. For hard sphere molecules, the viscosity 
is proportional to the square root of temperature and inversely proportional 
to the cross section (Ref. 4). Conversely, the effective cross section may be 
considered proportional to the square root of temperature and inversely pro- 
portional to viscosity. This is expressed as follows: 


(n-i) 


e - oo 

~a* 


e -oo 
T* , 


( 41 ) 
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where T e _ oc is an effective temperature corresponding to the relative velocity 
between colliding molecules, a* is the cross section at a reference temperature 
T* and n is the exponent in the power-law variation of viscosity with temper- 
ature. Based on molecular data in Ref. 4, a value of 43.0 x 10 cm for a 
at temperature T* = 300°K is assumed, and a value 0.8 for the exponent n. 

From atmospheric data in Ref. 3, = 41.9 x 10'“cn, 2 The effective temper- 

ature T can be defined in a number of ways, all of which should be based 
on a proportionality with the square of the relative velocity. We will define 
it as follows: 


T 

e -oo 


i jM 
3 R 

g 



oo 


(42) 


where M is the average molecular weight of the colliding molecules and R g 
is the universal gas constant. This expression is based on the mean-square 
velocity relationship to temperature for a Maxwellian gas. T e _ oo is thus 
found to be 50,200°K. Based on the above data, the cross section ratio 

o is found to be 0.22. 

c c 


Using the values which were obtained for Kn, S b and cr^^/cr^, the 

dimensionless factor defined by Eq. (12) is found to be 0.00665. This factor 

may be considered an inverse Knudsen number Kn^^ = RAg.*, characteristic 

of collisions between outflowing and ambient molecules. The characteristic 

mean-free-path X is seen to be much greater than the sphere radius, thus 
e -oo 

justifying the assumption of nearly free molecular flow. 


Using the above results and Eq. (15), the backscatter flux at * =0, due 
to interaction with the ambient atmosphere, is now found to be 0.003, or 0.3%, 
of the outgas rate. If the collision cross section had been considered constan 
and unaffected by the relative velocity, the cross section ratio “ ° ' 

CT =1.13, and the resulting backscatter flux would be 1.5% of the outgas rate. 

00 
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From Eqs . (35) and (40) it is seen that the fractional intensity of back- 
flow due to self-scattering is dependent on the outgas rate. First, assume 
an outgas rate based on an estimated 30 kg/day total weight loss rate for the 
Skylab spacecraft. Based on the 10 m radius for the equivalent sphere, this 
amounts to an outgas rate of 2.77 x 10 8 g/cm 2 /sec. Recall that this quantity 
includes waste dumps, attitude control reaction motor firings, etc., in addition 
to outgassing. This assumed value is actually about three orders of magnitude 
greater than that which could be expected for simple outgassing. 

From Eq. (35), the inverse Knudsen number Kn e _ e for self-scattering 
is found to be 1.39. The characteristic mean-free -path, therefore, is about 
0.67 of the sphere radius. This mean-f ree-path value is of marginally satis- 
factory magnitude for the first-collision assumptions used in the analysis. 

From Eq. (40), the backscatter flux due to self-scattering is found to be 0.30, 
or 30%, of the outgas rate. This value is quite large, and, in view of the 
lumping together of sucn contributions as waste dumps and reaction motor 
firings, the question arises as to the validity of the results. 

Waste dumps and reaction motor firings occur itermittently and at 
rather high instantaneous pressures compared to the average, assuming a 
continuous discharge. The flow regime, therefore, is actually much more 
nearly continuum than free molecular, and the mass iov. returning to the 
vehicle is more likely due to direct jet impingement than to mter molecular 
collisions. 

A typical outgas rate for spacecraft materials is 3.9 x 10 g/cm /sec 

for BBRC "one" paint, air dried (Ref. 5). For this outgas rate, and still 

assumin'. 1 M -24, Kn* is 2.4 x 10 3 . The mean-free -path is now about 400 

e-e 

times the sphere radius. This is well within the nearly free molecular flow 
regime. From Eq. (40) the backscatter flux due to self-scattering is found to 
be 0.0005, or 0.05%, of the outgas. rate. 
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Section 4 

CONCLUDING REMARKS 

The derivation of equations in this study is fairly rigorous for back- 
flow due to interaction with the ambient atmosphere. Some rather extreme 
simplifying assumptions were made, however, for the case of self-scattering. 
When outgassing alone is considered, the contribution to backflow due to self- 
scattering is found to be fairly small, about an order of magnitude small r 
than that due to interaction with the ambient atmosphere. The rather extreme 
simplifying assumptions for self- scatte ring, therefore, are found to be justified. 

When all forms of mass flow from the spacecraft, including waste dumps, 
reaction motor firings, etc., are averaged together into a single continuous out- 
flow rate, the density of this averaged outflow is sufficiently high that self- 
scattering is found !o far overshadow interactions with the ambient atmosphere. 
The simplifications made in the analysis of self-scattering, therefore, greatly 
reduce the accuracy of the calculated backflow . Also, as was pointed out in the 
preceding section, the primary mode of backflow onto the spacecraft from waste 
dumps and reaction motor firings is probably direct jet impingement rather than 
inte r molecular collisions. 

Considering outgassing alone, only a small fraction, i at most, of the 
outgas products return to the spacecraft as a result of inter molecular collisions. 
For some outgas products, even this small amount may significantly contaminate 
sensitive experiment packages. Considering the total outilow (waste dumps, re- 
action motor firings, etc.), a rather large fraction, about 30%, was estimated 
to return to the spacecraft. This is certainly a significant quantity; however 
questions of accuracy of the analysis, and even its applicability to waste dumps 
and jets, should be considered in evaluating the results. 
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FOREWORD 
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Section 1 

INTRODUCTION AND SUMMARY 

The expulsion of liquid wastes from an orbiting spacecraft results in 
the buildup of a cloud of ice particles which orbit about the central body 
(earth, moon, etc.) in much the same fashion as the spacecraft from which 
they were ejected. In fact, the trajectories of these particles may be not 
far removed from the orbit of the spacecraft, and some particles may even 
return to the spacecraft where the orbits intersect. There are two primary 
areas of concern about this cloud. One is the possibility of interference 
with optical experiments by light scattering and attenuation, and the other 
is the contamination of the spacecraft by those particles which eventually 
return to the spacecraft and stick or condense on it. 

To provide an indication of the rate of mass flow return to the space- 
craft, and the mass distribution of the particle cloud, a theoretical analysis 
was performed for the case of a spherical spacecraft in a circular orbit 
ejecting mass uniformly in all directions at a constant rate. The expulsion 
velocity was assumed identical for all particles. It was found that, neglect- 
ing aerodynamic drag, the rate of return to the Skylab spacecraft is of the 
order of the outgas rate of a typica’ spacecraft material. Aerodynamic drag 
was found to eliminate the possibility of return flow altogether for Skylab. 

The optical properties of the particle cloud were found to be not substantially 
more than that predicted for an isotropic source with t.ie >ame expulsion rate. 
Some intensification was found in directions along the sp? .recraft orbit. 
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Section 2 

DERIVATION OF EQUATIONS 


2.1 PARTICLE ORBITS 


The orbital motion of the ejected particles is determined in the follow- 
ing analysis by utilizing small perturbation theory as outlined in Ref. 1. The 
particle orbits are derived as perturbations from the circular orbit of the 
spacecraft. Spherical polar coordinates r, 0 and <£ are used, where r is the 
radial distance from the center of the spacecraft orbit, 9 is the angular dis- 
placement in the plane of the spacecraft orbit, and 0 is the angular displace- 
ment away from the spacecraft orbit plane. The coordinate system is illustrated 
in Fig. 1. Note that <t> is defined here somewhat differently than is customary 
with polar coordinates. 


The differential equations for the orbits are 


• • 2 i 

r - r cos </> 0 


cos <f> 


• • r * * 2 

(j) + 2 — <j> + 0 sin <f> cos</> 

r 


o 

II 

^ 1 

PI 

+ 

(i) 

m v 

lai . o 

(2) 

m v 

.1° i = 0 

m v 

(3) 


where /Lt = yM, with y the gravitational constant and M the mass ot the central 
body, F d is the aerodynamic drag for ;e, m is the particle mass and 

v = <i- 2 + r 2 cos 2 * i 2 + r 2 * 2 ) l/2 

is the instantaneous velocity. 
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l - Illustration of Coordinate System Used in Orbit Analysis 
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Taking small perturbations of Eqs. (1) and (2) about the circular spacecraft 
orbit yields: 


6r + 



m v 

o 


6r 



60 


0 


(4) 


2 — 6*r + r 6 *9 + - 0 

r o 111 

o 


(5) 


where the subscript zero refers to the spacecraft orbit, and 6 denotes a 
perturbation from the spacecraft orbit. Note that v q and r Q are constant 
for circular orbits. Equation (3) requires no small perturbation analysis 
since the angle <t> is itself a small perturbation from the spacecraft orbit 

<*c =■ 0) ’ 


The solution of the differential equations (Eqs. (3) to (5)) will require 
that initial values be known at the time of particle expulsion from the space- 
craft (t = 0). These initial values are listed as follows: 


v cos?, 
e 


ir t=0 = »■ = -T • 6r t=0 * v e sin & C0S "’ *t»0 


v sin£ sin/7 


( 6 ) 


where v is the expulsion velocity, and £ and 7] are spherical polar angles 
describing the direction of expulsion from the spacecraft. The angles £ and 
7 1 are defined in Fig. 2, and the initial values are shown in their relationship 

to t, and 7], 


Equation (5) may be integrated to yield: 

2r« r - 6r t=0 )+r o< 5 ' 9 - S ' 6 t=0> 


(7) 


With the appropriate initial values, Eq. (7) becomes: 


v • r 

Z 6r+r^69-v^ cos?, + t = 0 


( 8 ) 
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Relationship of Initial Values to Expulsion Direction 
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Equations (4) and (8) are now combined to eliminate 69: 


6r + 



m v 

o 


2 

v 

6r + — ^ 6r 
r 

o 


V V 

2 — — — cos £ - 2 


F D V o 
m r 


(9) 


Equation (9) is a rather simple linear nonhomogeneous equation with constant 
coefficients which can be solved with little difficulty. Before finding the 
solution, however, the order of magnitude of the terms in the equation will 
be estimated. The following is a list of order of magnitude estimates 
which we will use in this analysis. 

6r ~~ v T , 6r - v , 6r ~ v / T 
e o e e o 


F D ~ P co V 14’ m ~ P p d p 


( 10 ) 


where p^ is the atmospheric density at orbit altitude, is the mass density 
of the particle, d^ is the effective particle diameter and T q ~ Z i c r Q / v Q > t ^ ie 

spacecraft orbit period. 

For the 435 km altitude orbit corresponding to the Skylab spacecraft, 

12 3 6 

is approximately 3.86 x 10 kg/m , r Q is 6.8 x 10 m and v q is 

7640 m/sec. The effective particle diameter is estimated from a force bal- 
ance between surface tension and vapor pressure forces in a water droplet 
in vacuum. 


d 

P 


= 4 t /p v 


( 11 ) 
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where T is surface tension and p^ is vapor pressure. For water near the 
freezing temperature, T and p are approximately 75 dynes/cm and 6660 
dynes/cm (5mm Mg), respectively. The diameter d p is thus found to be 
about 0.450 mm. 

The expulsion velocity v is probably of the order of 10 m/sec, cor- 

^ / 2 

responding to a differential pressure of about 7000 kg/m 10 psi). The 

particle density p for water droplets is, of course, 10 kg/m . 

Based on the above results, the order of magnitude of the terms in 
Eq. (9) are determined and listed as follows: 

First term — v /T 
e o 

Second term ~ (Poo/PpH r 0 / d p )( v e / T 0 ) ~ 10 4 ( v e / T 0 ) 

Third term ~ v /T 
e o 

Fourth term ~ v /T 
(p /p )(r /d ) 

Fifth term - — — — (t/T )(v /T ) 

(v /v ) o e o 

v e o 

~ I0 '‘ <‘/T 0 > (T e /T o> 

The first, third and fourth terms are seen to be of about the same order cf 
magnitude. The second term is clearly negligible compared to the others. 

The fifth term becomes significant after a few orbit periods and progressively 
increases in magnitude thereafter. Actually, particle sizes can extend down to 
the micron range, thus increasing the relati e order of magnitude of the second 
term. Even for a half micron size particle, however, the factor is about 0.1, 
with the result that the second term may still be safely neglected. 
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Neglecting the second term, Eq. (9) is solved to yield: 

F T 

7—7 — r = (sin ; cos n + — -zrr 2 > sin < 2 * t ^ T o ) 

r o ( V\7 * m v e 


F T 

+ 2 cos C[l - cos (2?rt/T )] - 2 r °~ (t/T n ) 


( 12 ) 


m v 


or 


- (sin t, cos r) + a/~ ) sin (2 x t) 
+ 2 cos £ [l - cos (2 ?r ”t)] - 2 a t 


(13) 


where the tilde denotes normalization of 6r to the characteristic dimension 

r (v /v ), and t to the orbit period T . The aerodynamic drag parameter 
o e o o 

a is given by 


a 


F T 


o o 
m v 

e 


(14) 


Neglecting a, the variation in particle orbit radius from the space- 
craft orbit is illustrated in Fig. 3. For t, = 0, the particle orbit is tangent 
to the spacecraft orbit at the point corresponding to t = 0, 1, 2 . . .. For 
rj = 0 and E, / 0, the particle orbit crosses the spacecraft orbit at that 
point and some other point depending on the value of ?>. For 4 = if /2, the 

other point is at t = l/2, 3/2, 5/2, . . , For r) / 0, the projection of the 
particle orbit on the spacecraft orbit plane behaves in the same manner as 
for r) - 0. The effect of a is to continuously decrease the particle orbit 
radius. 
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Variation of Particle Orbit Radius from Spacecraft Orbit 
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Equation (12) may now be combined with Eq. (8) to yield 
v cos C Fpt v 

_ _ 3 _£ + 3 ■ — - 2 (sin£ cos r) + a /it ) sin(2 n t/T ) 

ao .. mr r ° 

o O o 


V 

f 4 — — cos C cos (2rt/T ) (15) 

r o 

o 


which is ;hen integrated and nondimensionalized to yield the final result for 
the variation of 69: 


L 


6 0 /( / v o ) - - 6 tt (cos C ) t + 37ra t 2 
- L (sinC, cost) + a/7r ) [ 1 - cos(2tt tf )] + 4 cos £ sin (27 r T ) 

1 

(16) ; 

4 

J 

| 

Again neglecting a, particles ejected ahead of the spacecraft ( t, < n /Z) are | 

seen to recede behind the spacecraft, and particles ejected behind the space- j 

craft ( £, v 7 r / 2) are seen to advance ahead of the spacecraft. Particles ejected 
normally to the spacecraft trajectory ( C, - tt /Z) and r) / + tt fZ are seen to 
advance and recede over the orbit period, so that 60 = 0 at the end of orbit 
periods. For C = tt /Z and rj = + jt/ 2, 60 remains constant at zero. 

The variation in <f> is found from the solution of the differential equation, 

Eq. (3). This equation is nonlinear in form, however, and cannot be readily 
solved. Before a solution is attempted, the order of magnitude of the .e-.ns 
in the equation will be estimated to determine if the equation can be realistic- 
ally simplified to a more tractable form. The following is a list of order of 
magnitude estimates to be used in estimating the order of magnitude of the 
various terms: 
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Also, the following near -equalitie s will be assumed: 


v = v , r = r 
o' o 


= r , 0 = — , sin^ = 0, cos^ =1 (<£ « 1) (18) 


Using these estimates, the order of magnitude of the terms in Eq. (3) are 
determined and listed as follows: 


First term ~ 


Second term 




r T \v I 2 
oo ' o T 

o 




jV \ v T / v \ i 
Third term ~ (^) 12 


Fourth term 


F^ v /V \ , 

— s S.~ o(_S.)_* 

mv o r o ' v o'T 2 

o 


The first and third terms are seen to be about the same order of magnitude. 

The seoond term differs by the factor v and the fourth term differs by 

C ° - 3 - 1 

a (v /v )„ Previous results show v /v ^ 10 and Of ^ 10 . The second 

v o *- o 

and tourth terms, therefore, will be neglected. Taking into account the near- 
equalities of Eqs. (18) and dropping the second and fourth terms, Eq. (3) be- 
comes: 


V + < v 0 / r Q ) ^ r - 0 


With the initial values in Eq. (6), the solution of Eq. (19) becomes: 
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$/(vV v o ) - sin 


E, sin rj sin (2 n t ) 


( 20 ) 


Equation (20) reveals that particles ejected from the spacecraft in a direc- 
tion out of the spacecraft orbit plane will swing back and forth across the 
orbit plane, crossing the plane at times corresponding to half and whole 
values of orbit periods (t = n/2, n = 1, 2, . . . ). 


2.2 RETURN MASS FLOW RATE 

Reviewing the results of the preceding section, neglecting a, the par- 
ticle position was found to coincide with the spacecraft position at the end 
of whole orbit periods when 4 = it/ 2 and at the end of half orbit periods when 
both Z, -■ it /Z and rj = +ir/Z. Therefore, particles ejected within two small 
solid angle regions about \ - it /Z and T) = it /Z will collide with a finite 
diameter spacecraft at the end of a naif orbit period, and particles ijected 
within a small range of £ about tt/2 (for all \) will collide with the space- 
craft at the end of a whole orbit period. In tne following paragraphs, the 
range of £ and T) for wuich collisions occur are found, and the resulting 
return mass flow rate is determined. 

The distance d of the particle from the spacecraft is given by (for 

d/r « 1): 
o 

1 ^2 

d/r Q = £(6r/r Q )^ + 6 9^ + \ (^1) 

By setting d equal to the spacecraft radius R (assumed spherical) at the times 
of closest approach (t = l/Z and 1), the limiting range of t, and n for particle 
collisions with the spacecraft can be determined. 


For t = l/2, we obtain the following expression: 


R = | (16 +9 it ) cos^ ^ + 24 r sin £ cos C cos rj 
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+ 16 sin 2 E, cos 2 17 - (b - 9 ;r 2 /l 6 - l 6 / 7 r 2 )ct 2 


+ [( 16 + 9 /r / 2 ) cos - (32 /tt - 67 t ) sin £ cost] J a 


|l/2 


( 22 ) 


For the Skylab spacecraft, an effective radius R of about 10 m may be 

assumed. Using the values of r , v and v given in the preceding section, 

o 0 o ^ 

the left side of Eq. (22) takes on the value 1.12 x 10 . This requires that, 

for Eq. (22) to be satisfied, the square root of each term under the radical 

-3 . -1 

be no greater than the order 10 . Since a is of order 10 for the Skylab 

spacecraft, this implies that no particle collisions can occur. Neglecting 
aerodynamic drag (a = 0), we find that cos C ^ cos rj ~ 10 f or C — 7T fl 
and rj == + 7 r /Z in accordance with our previous findings. Equation (22) may 
now be simplified by changing variables and using appropriate approximations 
as follows: 


(16 + 9tt 2 ) A 2 + 24 A. A + 16 A f = R 2 (23) 

t <s T] t 

where A r = tt / 2 - t, and A = + ir /Z - t] , Equation (23) is seen to be the 

^ f 

equation of very small ellipses about £ = 7 T / 2 and 77 = +_ tt /2 on a unit radius 
sphere. The solid angle enclosed by these ellipses is simply their area. 
From analytic geometry, the axis of the ellipse is seen to be rotated from the 
axis by an angle 9 given by 

0 = (1/2) tan " 1 (8/3 tt )= 20.17° (24) 


and the two semi -axes a and b are given by 


a 


R 

|(16+ 9 7 r 2 ) cos 2 9 + 24 it sin 0 cos 0 + 16 sin 2 9j 


1/2 
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= 0.0917R 


(25) 


and 



R 

24 7 T sin 9 cos 


9+ 16 cos 0 


1/2 


= 0.463 R ( 26 ) 

The area and, hence, the solid angle u is 

a) = 7r ab ( 2 ^) 

which, for the two ellipses, amounts to a fraction f of the total solid angle 
given by 


f = 2 7 r ab/ (47 t ) = ( l/ 2) ab 

= 0.02 12 R 2 ( 28 ) 

Using the value of R found in a preceding paragraph, this fraction is found 
to be 2.66 x 10’ 8 , which is the fraction of mass flowing from the spacecraft 
which will return a half o’-bit period later. 

For 7=1, the distance expression (Eq. (21)) becomes 

\ « / O 

R = | 36 7 T 2 cos 2 £ + (97T 2 + 4) a 2 - 367r 2 acos 
= | ( 6 7r cos^ - a ) 2 + ( 9 7r 2 + 3)a 2 | 
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As was the case for t = l/2, the expression cannot be satisfied for the aero- 
dynamic drag anticipated for the Skylab orbit, and, hence, no particle colli- 
sions are possible. Neglecting aerodynamic drag, we find that cos£ ~ 10 , 

or £ =7r/2, again in accordance with our previous findings. With the ap- 
propriate simplifications, Eq. (29) yields 

A ^ = +_ r/ 6 n (30) 

The solid angle enclosed by + A ^ is found by integrating as follows: 

A^ 

u = 4 7 r f cos A £ dA ^ sinA 4^ A ^ 

*T 

= (2/3) R (31) 

The fraction f of the total solid angle is 

f = R/ 6 7r (32) 

-5 

For the Skylab spacecraft, this fraction is found to be 5.94 x 10 , which is 

greater than the fraction for t = 1 / 2 by a factor of about 2000. Due to the 
relative insignificance of the return flow for t = l/2 compared to that for 
7=1, and since this fraction is included within the fraction for t = 1 anyway, 
we will consider only the fraction for t = 1 in our analyses for return mass 
flow rate. The return mass flow rate m can now be calculated from 

m = R fa/bn (33) 

r 

2.3 DENSITY FIELD ANALYSIS 

From the particle orbit equations derived in Section 2.1, the density 
field of particles ejected from the spacecraft over a period of time can be 
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determined. The density p at a point P and after a time t is found by carry- 
ing out the following integration: 


P P, t 


Lim 
6 V— * 0 


6 V 


fS 

o 4* 


rrx f y „ j dudt 

i,v,t 


(34) 


where m is the rate of mass expulsion from the spacecraft and f du is the 
fraction of this mass which is expelled during time interval dt within the 
solid angle du at £ and 77 , and which is contained within the volume element 
6 V at point P after a time t’ This expression includes in the density calcu- 
lation mass that is moving away from the spacecraft for the first time after 
having been ejected, and mass that may have been trapped within the particle 
cloud for some time. For two reasons, our primary concern in these calcu- 
lations is with the latter portion of mass. First, it is easy enough to calculate 
the density field for the source-like flow of mass moving away from the space- 
craft for the first time. Second, the singular nature of the density field near 
the source renders any numerical calculations meaningless where large volume 
elements are used. It will be shown later that, for light scattering purposes, 
the bulk of the density field of an isotropic source lies within a few source radii. 
For these reasons, we will not include in the numerical calculations that por- 
tion of mass in the source-like flow period of time. For a complete density 
field calculation, the source-flow density, which can be easily calculated by 
other means, can be added to the trapped particle density calculations. Since 
the point at which the source-like flow may be considered terminated is some- 
what arbitrary, we will choose a period of time in orbit of 0.3 of an orbit period 
(7=0.3) at which to begin our calculations. From Eq. (25), particles ejected 
out-of-plane reach the maximum out-of-plane displacement at t = 0.25. 

From Eqs. (13), (16) and (20) the maximum range of a particle in its orbit 
(neglecting orbit decay) is given by: 



(3 5) 
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60 


max 


+ 67T 


V 



v max 
o 


(36) 


* 


max 



(37) 


where t is the maximum number of orbit periods under consideration, 
max 

The 6r and <j) ranges are divided into 2N equal increments and the 60 range 

into 2N T increments to form volume elements of siz,e: 
max 



The number of 6 9 increments is made proportional to the number of orbits 

r because the orbit range is proportional to t . This keeps the volume 
max max 

element size the same regardless of the number of orbit periods to be con- 
sidered. 

Next, the total solid angle (4^ steradians) about the spacecraft is 
divided into nearly equal elements. This is accomplished by first dividing 
the C range from 0 to tt into 2M equal increments A C. Then, depending on 
the AC increment, the tj range from 0 to 2tt is divided into the appropriate 
number of increments such that the solid angle element Au enclosed by AC 
and At/ is about the same for each AC increment. This is illustrated in Fig. 4 
The solid angle element Au q between C = ® anc ^ £ = * s 

Au = 2^r ( 1 - cosAC ) ( 3< ?) 

o 

This conical-shaped element is not further divided into increments of rj but 
is taken as the basis for dividing the tj range for the other elements. The 
fraction ^ of the total solid angle included in Au q is 
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c^o = ~ir " i (1 ' cosA;) 


The solid angle AfT included between £ = i At, and t, = (i + 1 ) A£ is 

Afi. = 2ir jcos(iAU - cos [(i + 1 ) A&J j < 41 

This element is further divided into M,^ equal increments of i) such that the 

resulting elements Au. are about the same size as Aw q . Dividing AS^ by 

M,j and equating to Au q yields to the nearest integer: 
i 

M„. - Integer jcosjlAiI^plAil| ,41 

For most practical values of A£ , values are found to increase with i 

in the sequence 3, 5, 7, . . . for the first few values of i. The fraction <^7 of 
the total solid angle included in Au^ is 


04 i = ~47 


cos(iA^) - cosL(i + 1 ) A( 
2 M„ 


Recall that these values of should be close to 

Finally, using Eqs. (13), ( 1 6) and (20), an amount of mass equal to 
mAt & . is considered to be ejected from the spacecraft at an instant of time, 
and a determination is made as to the volume element in which the mass of 
particles is contained at later times. This is equivalent to solving Eq. (34) 
in the following fashion: 


m At 


N t 2M 


'i j k'N t ‘ AV 




n t =l i=l j=l 
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r i > 

where i and j refer to the £, and 17 increments, respectively; i , j and k 
refer to the volume element corresponding to the 5r, 60 and <j> increments, 
respectively; n fc refers to the time step; and ^ j i'j Vn t = 1 or 0 depending 

on whether the mass of particles from the i j solid angle element is or is not 
respectively, in the i' j'k volume element after a time interval n fc At + 0.3 T Q . 
(Recall that the time steps begin at t = 0.3.) Including the expression for AV 
from Eq. (38), and normalizing At to the orbit period T q , yields; 


^Vjk 


24ff 


m T At 

o 

(v /v ) r 
e / o o 

N 


EEE 


(45) 


For convenience, we normalize the local density to an average density p for 
mass ejected during a single orbit period and over the maximum range of 
particles during a single orbit period. 


P 


m T 


o 




(46) 


The final result for the normalized density expression is now 

Pi'jVN iii, 2M “j! 

^jVN t * — W~ - 8n3aT E Z S 

f n t = l i = l j = l 


(47) 


A computation procedure based on Eq. (47) was programmed for the 
NASA-MSFC Univac 1108 digital computer system. A listing of this program 
is given in the Appendix. 


2.4 COLUMN DENSITIES 

The parameter of interest in light scattering by the p 9 r k i\ cloud is the 
column density along a line of sight. This paramete. is determined as follows: 
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p(s) ds 


(48) 


where N is the column density and the integration is taken along the line of 
sight r. In the present study, the line of sight will originate at the spacecraft 
and extend out in some direction of interest beyond the outer extent of the 
particle cloud. We will compare the calculated column densitie s for the 
particle cloud with a theoretical column density for an isotropic spherical 
source of radius R. This theoretical column density is given by 


N 



R 


m 


dr 


4 7T r 


m 

4 7T R v 

e 


(49) 


where r is the radial distance out from the center of the source. Note from 
Eq. (49) that the value of the integrand continuously decreases and approaches 
zero as r increases. There should be some characteristic distance R' out 
from the source, therefore, within which is contained the bulk of the integral 
value. By integrating out to the limit R' and requiring the integral value to be 
some fraction, say 0.9, of the total integral, we determine that 90% of the 
column density value for an isotropic source is contained within a distance 

of 10 source radii. 

Normalizing the column density N c in Eq. (48) to the isotropic source 
column density N' c in Eq. (49), and normalizing the source radius R and the 
distance s along the line of sight to the characteristic length ( v e / v 0 ) r 0 for 
the parficle cloud, we change Eq. (48) to the following form. 


N c =£R / p(s)d s 
T 

where the tilde denotes normalized values. Equation (50) may be integrated 
numerically along a line of sight originating at the spacecraft location and 
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extending in any direction. The computer program listed in the appendix 
performs the integration in six different directions, i.e., in the positive and 
negative r, 6 and <f> directions away from the spacecraft. The integration 
procedure uses the trapezoidal formula with the step sizes identical to that 
used in forming volume elements for the density field analysis. 
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Section 3 
RESULTS 

3. i RETURN MASS FLOW RATE 

The return mass flow rate to the Skylab spacecraft was estimated based 
on the parameters determined in Section 2.2 and on an estimated total weight 
loss rate of 30 kg/ciay. Neglecting aerodynamic drag, the return mass flow 
rate is determined from Eq. (33) to be approximately 2x 10 kg/sec. (Recall 
that the predicted effect of aerodynamic drag is to prevent mass return to the 
Skylab.) Based on the 10 m effective radius, this is equivalent to a flux of 
about 0.7 x 10 -11 g/cm 2 /sec. A typical outgas rate for spacecraft materials 
is 4 x 10' U g/cm 2 /sec for BBRC "one" paint, air-dried (Ref. 2). The pre- 
dicted return mass flow rate is seen to be of the same order of magnitude as 
this outgas rate. 

3.2 PARTICLE CLOUD DENSITY FIELD 

The computer program based on Eq. (47) was used to calculate the norm- 
alized particle density field about an orbiting spacecraft for durations in orbit 
up to 10 orbit periods. The dimensionless character of the results makes them 
applicable to arbitrary conditions of mass expulsion rate and orbital parameters. j 

Calculations were made for zero aerodynamic drag and for drag conditions cor- , 

responding to that anticipated for water droplets expelled from the Skylab space- 
craft. A particle orbit decay parameter ct value of 0.461 was used where drag 
was considered. Equation (14) was used to calculate a, with the drag force F D 
calculated using a drag force coefficient of 2.2, and with the particle diameter, 
atmospheric density, etc., the same as those used in Section 2. 

Results of the density field calculations are shown in Figs. 5 through 12 
lor a time in orbit up to 10 orbit periods. The position coordinates are appro- 
priately normalized for arbitrary applicability to various expulsion velocities 
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and spacecraft orbit parameters. The mesh parameters used in these calcu- 
lations were M = 10 and N = 5, and the time step A? was 0.01. Recall that 
these calculations do not include particles in the source -flow phase of their 
trajectories (moving away from the spacecraft for the first time after having 
been ejected). The density field for an isotropic source is shown separately 
on the figures for comparison with the calculated results. The total density 
in the particle cloud is the sum of the two densities. 


The data points in Figs. 5 through 8 did not show very much scatter, 
and so were simply connected by straight lines to display the overall trend. 

In Figs. 9 and 10, however, considerable scatter is present, and the trends 
are indicated by coarsly fairing the data using straight-line segments. This 
scatter in the data is evidently the result of numerical error stemming from 
the finite coarseness of the mesh and the finite time steps. By experimenting 
with the mesh and time step sizes, it was found that the scatter could be re- 
duced. The computational costs increased, however, to the point that the 
increased accuracy was not economically justified. And, even with the data 
scatter, the overall trends seem quite clear and reasonable. 


Figures 5 and 6 show that the density distribution in the radial direction 
(6 ) maximizes and reaches an equilibrium after about six orbit periods for both 
the zero drag and finite drag cases. Aerodynamic drag is shown in Fig. 6 to 
shift the distribution toward the central body (-6r direction). 


The density distribution in the out-of-plane (<M direction is shown in 
Figs. 7 and 8. Symmetry in the (f> direction is shown to be unaffected by aerody- 
namic drag. The density level, however, is seen to be shifted down and smoothed 
out somewhat by the effects of drag. The downward shift at the spacecraft 
location is due to the overall shift of the cloud in the negative radial direction 
(Fig. 6). The slight peaking toward the edge shown in Fig. 7 is probably a com- 
putational anomaly caused by the finite grid mesh and time step. 
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The zero drag density distribution along the spacecraft orbit (6 0 
direction) is shown in Fig. 9 to maximize and reach an equilibrium at the 
spacecraft location; however, the extent of the distribution continues to ex- 
pand in the positive and negative 5 6 directions. These observations are in 
accord with the behavior of the particle orbit equations. The effect of drag, 
as seen in Fig. 10, is to cause an equilibrium to be reached in the expansion 
of the cloud along the spacecraft orbit. Actually, the cloud will continue to 
expand, but it will simply spiral inward toward the central body and away 
from the spacecraft orbit. 

The particle cloud shape was found to be symmetrical in the out-of-plane, 
<f > , direction at all points ahead of and behind the spacecraft on the spacecraft 
orbit. As shown in Fig. 11 for zero drag, however, the shape is skewed with 
respect to the radial direction. Ahead of the spacecraft, the cloud shifts in- 
ward in the negative radi?l direction, while behind the spacecraft, the cloud 
shifts outward. Based on Eqs. (35) through (37), the particle cloud about the 
Skylab spacecraft (for zero drag) is found to extend, on both sides of the space- 
craft, 35.6 km in the radial direction, 8.9 km in the out-of-plane direction and, 
after 10 orbits, 1680km along the spacecraft orbit. In the absence of drag, the 
particle cloud would extend completely around the orbit and begin to overlap at 
the spacecraft location after 255 orbit periods. 

A portion of the cloud shape for the case of aerodynamic drag is shown 
ir. Fig. 12. Calculations were not made beyond the lower extent shown. 

3.3 COLUMN DENSITIES 

Calculated values of column densities, normalized to the theoretical 
isotropic source value, are shown in Figs. 13 and 14. A value of R = 1.12 x 
10’ 3 is used which is based on the Skylab spacecraft. These results are for 
lines of sight originating at the spacecraft location and extending in both direc- 
tions along the r, 0 and <f> axes. The curvature of the orbit is not taken into 


3-11 


LOCKHEEt. HUNTSVILLE RESEARCH & ENGINEERING CENTER 




LMSC-HREC D306222 



r| = 180 



= 0 ° 

& 270 


I J 

5 . . ~ 10 
Normalized Time in Orbit, t 

Time Variation of Column Density Along Three Axes 

(Aerodynamic Drag Included) 




LMSC-HREC D306222 


account in these calculations. Also, the calculated column density values 
are based on the previously discussed cloud density calculations which do not 
include particles in the source flow portion of their trajectories. 

Figure 13 shows that the zero drag column densities in the radial and 
out-of-plane directions (r and 0) reach an equilibrium after about five orbits. 
The equilibrium values are seen to be of the order of 10 of the theoretical 
isot-opic source value. The column density in the direction along the orbit 
(0) continues to increase as the cloud expands along the orbit. The curvature 
of the orbit, however, which is not taken into account in these calculations, 
will eventually bring about an equilibrium value in this direction also. For 
the Skylab spacecraft, the maximum line of sight distance, based^on the cloud 
dimensions, is found to be about 2000 km. Using the maximum p value of 
27.5 in Figs. 5,7 and 9 as the average along the linejrf sight, Eq. (50) yields 
an estimated maximum normalized column density N £ of about 0.9. The actual 
value is certainly somewhat less than this estimation. 

When drag is considered, the calculated column densities in all direc- 
tions reach an equilibrium after about five orbits. These values are all seen 
to be of the order of 10" 2 to 10' 3 of the theoretical isotropic source value. 

The results of the column density calculations show that the buildup of 
the particle cloud in orbit does not significantly affect column densities and, 
hence, light scattering and attentuation effects about the spacecraft. As pointed 
out in Section 2.4, the bulk of the column density for an isotropic source is con- 
centrated within a few source radii, and the trapped particle cloud adds very 
little to that. For the case of zero drag, the trapped particle cloud effects are 
seen only within a narrow range in directions along the spacecraft orbit. 
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| Section 4 

CONCLUSIONS AND RECOMMENDATIONS 

^ The results of this study effort show that a particle cloud does indeed 

\ build up about an orbiting spacecraft due to waste dumping, and a small fraction 

1 of this mass may return to the spacecraft after times corresponding to one half 

■ and one orbit period. The estimated return rate to the Skylab spacecraft, neg- 

I lecting aerodynamic drag, was found to be of the same order as the outgas rate 

of a typical spacecraft material . When drag was considered, however, the 
| particle orbits were found to decay away from the spacecraft orbit such that 

no collisions with the spacecraft were possible. 

^ The trapped particle cloud was found to contribute very little to 

j light scattering and attenuation effects, as predicted by calculated column 

' densities, when compared with the source Dow region near the spacecraft, 

t The bulk of the column density along a line of sight originating at the space- 

I craft is concentrated within a few spacecraft radii. Some intensification of 

these effects, however, is found in directions c.long the spacecraft orbit. 

The analysis reported herein utilized certain simplifying assumptions, 
j such as a uniform distribution of particle sizes and expulsion velocities, a 

constant and uniform expulsion rate and spherical symmetry in the direction 
j of expulsion. Moreover, the particles were not considered to evaporate or 

otherwise be altered over a period of time. Future effort in this area should 
. include removing some of the restraints in these assumptions and possibly 

J. considering the long-term interaction and thermal equilibrium of these par- 

ticles with the ambient atmosphere. 
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Appendix 

j 

The computer program used in the particle cloud analysis is listed in 
j the following pages. This program was developed for use on the NASA -MS FC 

‘ Univac 1108 system. The input and output variables are described as follows: 

Input Variables 


I 

I 

I 


M 

M in Section 2.3 

N 

N in Section 2.3 

A 

a in Section 2 . 1 

TMAX 

T in Section 2.3 

max 

DELT 

At in Section 2.3 

TWRT 

instructs program to write output 
variables at t = u TWRT, n = 1, 2 

R 

R in Section 2 .2 


I 

l 

l 

I 

I 

I 

I 


Output Variable s 


F 

in Section 2.3 

FO 

in Section 2.3 

I 

i in Section 2.3 

IP, JP, KP 

i' ^ in Section 

DEN 

p in Section 2.3 

ME 

Mjj in Section 2.3 

T 

7 in Section 2.1 

GOLDEN 

N in Section 2.4 

c 

ETA 

if in Sect* on 2 .1 

ZETA 

£ in Section 2.1 
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PROGRAM LISTING 

otmfmfion f( ioo) *mf< ion> 

COMMON T«TMAX»ST .ct.fxpt. a.fz.cz.se.cf.pi .dr.dth.phi 
*« rEN( 10. 100* 10 ) « N * NP * F I *DFLT*N?N , N?NP . OLDNP « DEL *OELJ 

* « ) *PM ax » PPM IM.THTMAX.THTMIN* PH ! M AX « PH I M I N 

r?t An ( f * io«FNn=R°p)M,N * a > tmax *oflt • twRt . r 

WRlT'-( 6 »? 0 )M,N, A.TMAX.nFLT.TWRT.R 
in FOPf'ftT ( I = , JF • F 1 o. A « 3F 1 n . 3 .F10.4) 

20 FORM aT ( //^X?HM= I 3 , ?X?HN= J 3 . ?X?HA=F 10.4. 2X5HTMAX=F7. 2 » 2XFHDELT 

* F5.7.2XFH rwRT = er F.?* 2 X?HR = Fl 0.4 ) 

PI = 7, 1 4 lORPftF 
TCK=TWRT-O.F*nELT 

NP- I F I X ( TMAX*FLOAT ( N ) ) 

N°N=?*N- 1 
N2NP=?#NP- 1 
mo = 2 *m 

O' ,DNP=n.O*FLOAT ( N**3) *DFLT 
po-L7 = PI /FLOA~ ( ?*(M+1 ) ) 

D"L= 1 . O/^LOAT ( N ) 
n<TLJ=l • O /FLOAT ( NP ) 

R *M A y= ( FLOAT C N ) - <> . ft ) *OFL 
P1 »MTn = -RPMAx 

THTM AX= ( FLOAT ( NP ’ -0.5 ) *DFL J 
TMtM 1 N=*-THT M AX 

pi iimax=Rpmax 

PI 1 1 M I N~— PH 1 MAX 
CIV I = C0F (OFL7 ) 

Fo=o,^* ( i . o-cof 1 ) 

WR1 TP ( ft* 37 )FO 

33 FORMAT! //?X7HFO=F1 7.4 ) 

no 7n ! = 1 *M 

COFfol =COF (FLOAT ( 1+1 )*PPL7) 

XI =Cnc 1 -C0F1P1 

X7=X1/<2.0*FO) 

MF( I ) = I FIX' <?) 

Mt T = T F T X ( X2 + 0 . c ) 

Ir (MrT.GT.MF! I ) > MF ( I ) =MFT 
F ( I )=X 1 /FLOAT (?*MF( I ) ) 

30 COSI=COSIP1 
WR I Tr ( ft « 34 ) 

34 F )RMaT( // 4X1HI * 3X2HMF *ftXl HF// ) 

UlR I TF < ft . 7ft > ( 1 ♦ MF < I ) » F ( I) « ! = 1 * M ) 

3r FORMAT! 21ft.Fl 2.4 ) 

Trn . 3 

4 n 1 = 1* N2N 
00 40 J=1«N2NP 
no 40 K = 1 ,N2N 

40 OEN ( I , J»K) =0.0 

100 TP-?. 0 *Pl*r 

F'.*=FtN( TP) 

CT=COF( TP) 

FXPT =FXP ( — A*T ) 

F!=Fo 
F7=0.0 
C-T=l ,n 
F*“ = 0,0 

cr= 1 ,0 
CLL POf 
C*.LL OFNC 
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no ?oo 1= l «mp 

7= ( FLOAT ( ! ) +0 » 5 ) #DFL7 
57=5 fN< 7 ) 

'-7=C04< 7) 

!F( 1 .OT.MIOO TO 150 
MFP=MF< I ) 

F ! =F | ) 

OO TO 160 
150 MtPP = MP-I+l 
MFP=MF(MFPP) 
fi=f(mfpp) 

160 DFLF=P.0*PI/FL0AT(MFP) 

DO POO J=1.MFP 
F = ( FL O AT ( J 1 -0 . = ) *DFLF 
5F= 5 1 N < F 1 
CF=C05(F) 

C*,LL P05 
POO CALL DFN5 
FT=Fo 
57-0.0 
C2=— 1 .0 
5r;= o.o 
OF = 1 .O 
CALL POS 
call DFN5 

IF( T.LT.TCK )G0 TO 400 
TOK=TO<+TWPT 
'>f r, ITF(6»P60)T 
?50 FORMAT ( //10X?HT=FR»?) 

DO TOO JP= 1 . N?NF 
WT I TF ( 6 • ?5 1 ) JP 
?51 F.'RMAT ( /5XTHJP= T 0/ ) 

D<» Too <P=1,N?N 

TOO WW I TF ( 6 * TO 1) ( DFN ( I P * JP . KP ) . I P= 1 * N2N ) 

TO 1 FORMAT! 1 OF! 2.4 ) 

CtiPl*R / 24 « O/FLOAT < N ) 

C ! =C*6.0*P I 

COlDFN=C1. /?.0*DFN< N.NP.N ) 

N! 'PrMP+l 

no 600 J=NPP.N?NP 
600 COLDFN=COLDFN+C1 *DFN(N. J,N) 

WRI TF( 6.60P 1 COLDFN 

60? FORMAT ( //6y6H7FTA = 0, ?X7HC0LDFN=F 1 P.4 ) 

CCLDfn=C 1/P.o*DFN( N.NP.N) 

NPP=NP- 1 
DO 610 J=1,NPP 

610 CULDFN=COLDFN+Cl *DFN(N, J.N) 

Wo I TF (6.61? ) COLDFN 

61? FCRMAT(//6X1?H7FTA=1R0 DFG . ?X7HC0LDEN=F 1 2 .4 ) 

Cl =4,0*0 

C0LDFN=C1/?.0*DFN< N.NP.N) 

N(>P = m+ t 

DO 6^0 J=NPP.N?N 
6?o COLDFN=COLDFN+Cl *DFN( l.NP.N) 

WP l TF < 6 « 6PP ) COLDFN 

62P F^>RMaT(//FX 17HZFTA = O0 OFG » FT A = 0 » ?X7HC0LDFN=F 12.4) 
COLDFN=Cl /?«0*DFN< N.NP.N) 

NPP=M- 1 
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DO 600 I=1,NPP 

630 COLDFN=COLDFN+C1 *DFN( I ,NP.N) 

W"' I TP ( 6 » 63? ; COLDFN 

63? FORMATf //5X?3H7FTA=R0 DFG * FT As 1 flO Dr G. ?X7HC0LDFN=F 1 2.4 ) 
C.'LDFN=C/2.0*DFN (N*NP *N) 

N|>P=M+1 

DO 640 K=NPP.N?N 
640 CuLDfn=COLDFN+C*DFN(N,NP,K) 

Wl M TP ( 6 . 642 ) COLDFN 

642 F"RMAT( //5X?2HZETA=90 DFG*FTA=90 DFG ♦ 2X7HC0LDEN=E 1 2 . 4 ) 
COLD*^N=C/2 • 0*DFN ( N « NP ,N ) 

KH >P = N-1 

DO 6«0 K=1 i NPP 

650 COLDfn=COLDFN+C*DEN(N«NP«K) 

WR I Tf ( 6 , 652 ) COLDFN 

652 FORMAT ( //5X23HZETA =90 DFG* FT A = ?7C DFG. 2X7HC0LDFN=F 1 2.4 ) 

400 T = T +DFLT 

IF ( T.LT.TMAX) GO TO 1 OO 
QQQ GTOP 
END 

F'JBRoutINF POF 

COMMON T.TMAX.ST.CT.FXPT. A.SZ.CZ.SF.CF.PI .DR.DTH.PHI 

* « DEN ( 10* 100. 10) *N*NP«FI . DELT * N2N * N2NP ,DLDNP . DEL .DEL J 
*«'?PMAX.RPMIN.THTMAX.THTM!N«PHIM4X*PHIMIN 

DR=0.?F-*( ( S7*CE+A/P! ) *ST+?.0*C7*( 1 . 0-CT ) -?• 0*A*T ) 

D~H=-CZ*T/TMAX+0 ,f*A*T**?/TMAX- 1 .0/( 3.0*PI *TMAX)* ( ( SZ*CF+A/P I )* 

* f 1-CT)-?.0*C7*5T) 

r>| |I =C7#GF*5T 

RETURN 

c NO 

FliRPoUTlNF DFN5 

COMMON T .TMflX * ST ,CT *FXPT . A « SZ . CZ« SF * CE »PI *DR»DTH«PHI 
*.DEN( 10* 10O, 10) «N«NP»FI *DFLT,N2N*N2NP.DLDNP.DEL*DFLJ 
*,RPmAX.RPMIN,THTMAX.THTMIN,PH!MAX.PHIMIN 

LOGICAL LI .L2.L3 

LI =DD. GT.RPMAX.OR.DR.lt. PPM IN 
L' 5 = DTH.GT.THTMAX.0R.DTH.LT.THTMTN 
L 3 = P h I .GT. OH I MAX. OP. PH I .LT.PH1MIN 
! r (Ll . OR . L 2 • OR • L 7 ) GO TO 1 0O 
T CKP = TM A Y— . 5*DFL T 
I*"*! T .GT.TC<P)G0 TO 5 
I *= IFIX( T) + l 
n!t1 sNP-N* IT+ 1 
N!T?=NP+N*IT-1 
Gl> TO 
5 N*T1=1 

N I T2 = N?NP 
fl CONTINUE 

D< • 1 O I = l « N2N 
RP= ( FLOAT ( t — N ) +0 ,F ) *DFL 
IF(DO.GT.RP) GO TO 10 
tP= I 

GO TO 1 5 
in CONTI NUF 
IF oo ?o J=N I T 1 .NIT? 

THP= (FLOAT ( J-NP ) 4-0 . 5 ) *DFI — I 
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I K ( DTH • GT • THP ) GO TO 20 
JU= J 

GO TO 25 
20 CONTINUE 
25 DO 30 K = 1 »N2N 

PH ! P= ( FLOA T ( K-N ) 4 O . 5 ) *DFL 
I r* ( PH ! • GT • PH I P ) GO TO 30 
KP = K 

GO TO 35 
30 CONTINUE 

35 DEN ( IP, JP,KP) =DEN< I P , JP , KP ) +DLDNP*F I 
100 CONTINUE 
R c TtJON 
END 
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